libraries

library(tidyverse)
library(brms)
library(bayesplot)
library(here)
library(glue)
library(scales)

theme_set(theme_light())

source(glue("{params$common_dir_str}/brms_model.R"))
source(glue("{params$model_dir_str}/model_prior.R"))

data

obs_only <- 
  read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
  mutate(subj = as_factor(subj),
         obs_degree = error,
         error = obs_degree * (pi/180))
## Parsed with column specification:
## cols(
##   subj = col_double(),
##   subj_index = col_double(),
##   stimulation = col_double(),
##   error = col_double()
## )

full model

print(bprior_full)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.5)     b   intercept            circSD            
## 2   normal(0, 0.5)     b stimulation            circSD            
## 3   normal(0, 0.5)    sd   Intercept  subj      circSD            
## 4   normal(0, 0.5)    sd stimulation  subj      circSD            
## 5   normal(0, 1.5)     b   intercept             theta            
## 6   normal(0, 1.5)     b stimulation             theta            
## 7     normal(0, 1)    sd   Intercept  subj       theta            
## 8     normal(0, 1)    sd stimulation  subj       theta
iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

fit_full <- brm(bform_full, obs_only, prior = bprior_full, 
                family = vm_uniform_mix, stanvars = stanvars,
                warmup = warmup, iter = iter, cores = cores, chains = chains, 
                control = list(adapt_delta = 0.99), inits = 0, 
                file =  glue("{params$save_dir_str}/fit_full"))

fit check

divergences

model_fit <- fit_full

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on circSD parameter

print(bprior_DcircSD_null)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.5)     b   intercept            circSD            
## 2   normal(0, 0.5)    sd   Intercept  subj      circSD            
## 3   normal(0, 0.5)    sd stimulation  subj      circSD            
## 4   normal(0, 1.5)     b   intercept             theta            
## 5   normal(0, 1.5)     b stimulation             theta            
## 6     normal(0, 1)    sd   Intercept  subj       theta            
## 7     normal(0, 1)    sd stimulation  subj       theta
fit_DcircSD_null <- brm(bform_DcircSD_null, obs_only, prior = bprior_DcircSD_null, 
                        family = vm_uniform_mix, stanvars = stanvars,
                        warmup = warmup, iter = iter, cores = cores, chains = chains, 
                        control = list(adapt_delta = 0.99), inits = 0, 
                        file =  glue("{params$save_dir_str}/fit_DcircSD_null"))

fit check

divergences

model_fit <- fit_DcircSD_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on pMem parameter

print(bprior_DpMem_null)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.5)     b   intercept            circSD            
## 2   normal(0, 0.5)     b stimulation            circSD            
## 3   normal(0, 0.5)    sd   Intercept  subj      circSD            
## 4   normal(0, 0.5)    sd stimulation  subj      circSD            
## 5   normal(0, 1.5)     b   intercept             theta            
## 6     normal(0, 1)    sd   Intercept  subj       theta            
## 7     normal(0, 1)    sd stimulation  subj       theta
fit_DpMem_null <- brm(bform_DpMem_null, obs_only, prior = bprior_DpMem_null, 
                      family = vm_uniform_mix, stanvars = stanvars,
                      warmup = warmup, iter = iter, cores = cores, chains = chains, 
                      control = list(adapt_delta = 0.99), inits = 0, 
                      file = glue("{params$save_dir_str}/fit_DpMem_null"))

fit check

divergences

model_fit <- fit_DpMem_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on pMem or circSD parameter

bprior_DcircSD_DpMem_null
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.5)     b   intercept            circSD            
## 2   normal(0, 0.5)    sd   Intercept  subj      circSD            
## 3   normal(0, 0.5)    sd stimulation  subj      circSD            
## 4   normal(0, 1.5)     b   intercept             theta            
## 5     normal(0, 1)    sd   Intercept  subj       theta            
## 6     normal(0, 1)    sd stimulation  subj       theta
fit_DcircSD_DpMem_null <- 
              brm(bform_DcircSD_DpMem_null, obs_only, prior = bprior_DcircSD_DpMem_null, 
                  family = vm_uniform_mix, stanvars = stanvars,
                  warmup = warmup, iter = iter, cores = cores, chains = chains, 
                  control = list(adapt_delta = 0.99), inits = 0, 
                  file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))

fit check

divergences

model_fit <- fit_DcircSD_DpMem_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

loo compare

#need to do this only once?
expose_functions(fit_full, vectorize = TRUE)

fit_full <- add_waic(fit_full, file = glue("{params$save_dir_str}/fit_full"))
fit_DcircSD_null <- add_waic(fit_DcircSD_null, file = glue("{params$save_dir_str}/fit_DcircSD_null"))
fit_DpMem_null <- add_waic(fit_DpMem_null, file = glue("{params$save_dir_str}/fit_DpMem_null"))
fit_DcircSD_DpMem_null <- add_waic(fit_DcircSD_DpMem_null, file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))

print("fit_full")
## [1] "fit_full"
fit_full$waic
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -804.9 15.9
## p_waic         7.2  0.2
## waic        1609.8 31.8
print("fit_DcircSD_null")
## [1] "fit_DcircSD_null"
fit_DcircSD_null$waic
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -805.3 15.9
## p_waic         7.2  0.2
## waic        1610.5 31.7
print("fit_DpMem_null")
## [1] "fit_DpMem_null"
fit_DpMem_null$waic
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -804.1 15.9
## p_waic         6.5  0.2
## waic        1608.3 31.8
print("fit_DcircSD_DpMem_null")
## [1] "fit_DcircSD_DpMem_null"
fit_DcircSD_DpMem_null$waic
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -804.6 15.9
## p_waic         6.7  0.2
## waic        1609.3 31.7
loo_compare(fit_full, fit_DcircSD_null, fit_DpMem_null, fit_DcircSD_DpMem_null, criterion = "waic")
##                        elpd_diff se_diff
## fit_DpMem_null          0.0       0.0   
## fit_DcircSD_DpMem_null -0.5       0.6   
## fit_full               -0.8       0.2   
## fit_DcircSD_null       -1.1       0.7

*Notes

All models were fit well.

WAIC diffs were not very large and had non-neglible SE. Also inconclusive.

waic diff ordering:

fit_DpMem_null

fit_DcircSD_DpMem_null

fit_full

fit_DcircSD_null

fit_DpMem_null predicted best since it omits the difference parameter on pMem, for which there seems little evidence of an effect.

fit_DcircSD_null fit worst since it omits the parameter that seems to be actually showing a non-zero value.

fit_full and fit_DcircSD_DpMem_null are in the middle.